Constraint-based, Single-point Approximate Kinetic Energy Functionals 
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We present a substantial extension of our constraint-based approach for development of orbital-free 
(OF) kinetic-energy (KE) density functionals intended for the calculation of quantum-mechanical 
forces in multi-scale molecular dynamics simulations. Suitability for realistic system simulations 
requires that the OF-KE functional yield accurate forces on the nuclei yet be relatively simple. 
We therefore require that the functionals be based on DFT constraints, local, dependent upon 
a small number of parameters fitted to a training set of limited size, and applicable beyond the 
scope of the training set. Our previous "modified conjoint" generalized-gradient-type functionals 
were constrained to producing a positive-definite Pauli potential. Though distinctly better than 
several published GGA-type functionals in that they gave semi-quantitative agreement with Born- 
Oppenheimer forces from full Kohn-Sham results, those modified conjoint functionals suffer from 
unphysical singularities at the nuclei. Here we show how to remove such singularities by introducing 
higher-order density derivatives. We give a simple illustration of such a functional used for the 
dissociation energy as a function of bond length for selected molecules. 

PACS numbers: 71.15.Mb, 31.15.xv, 31.15.E- 
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I. INTRODUCTION 

Simulation of the structure and properties of compli- 
cated materials is a demanding task, particularly away 
from equilibrium, for example, in the simultaneous pres- 
ence of solvents and mechanical stress. Though the 
present results are not limited to it, our motivating prob- 
lem has been tensile fracture of silica in the presence of 
water. 

For such problems, quantum mechanical treatment of 
the reactive zone is essential at least at the level of re- 
alistic Born-Oppenheimer (B-0) forces to drive an oth- 
erwise classical molecular dynamics (MD) or molecular 
mechanics (MM) calculation. Computational cost then 
leads to a nested-region strategy. Internuclear forces in 
the reactive zone are obtained from an explicitly quan- 
tum mechanical treatment. Forces between other nuclei 
are calculated from classical potentials. Such partitioning 
is called multi-scale simulation in the computational ma- 
terials community and QM/MD (or QM/MM) method- 
ology in computational molecular biology. 

The QM calculation is the computationally rate- 
limiting step. QM approximations good enough yet com- 
putationally fast on the scale of the MD algorithms are 
therefore critical. Both the inherent form and the grow- 
ing dominance of density functional theory (DFT) for de- 
scribing molecular, bio- molecular, and materials systems 
make it a reasonable candidate QM. Despite advances 
in pseudopotentials and order-A^ approximations, how- 
ever, solution of the DFT Kohn-Sham (KS) problem is 
too slow computationally to be fully satisfactory. An al- 
ternative, Car-Parrinello dynamics^, does not guarantee 
that the motion is restricted to the B-0 energy surface. 

Thus there is continuing need for methods which 



yield essentially full DFT accuracy at significantly lower 
computational cost. In response, we and co-workers 
proposed and demonstrated a Graded Sequence of 
Approximations^ scheme. Its essence is use of a sim- 
ple, classical potential for the majority of MD steps, 
with periodic correction by calibration to forces obtained 
from more accurate but slower methods. An example 
Graded Sequence of Approximations would be: (1) classi- 
cal potential; (2) simple reactive (charge re-distribution) 
potcntial^iii^; (3) Orbital-free (OF) DFT, the subject 
of this paper; (4) Quasi-spin density DFT^ (a way to 
approximate spin-dependent effects at the cost of non- 
spin-polarized KS-DFT); (5) FuU spin-polarized DFT 
(the level of refinement ultimately required for bond- 
breaking) . 

In this hierarchy, a large gap in computational cost sep- 
arates reactive potentials and quasi-spin density DFT. 
Since the cost of conventional KS calculations comes 
from solving for the KS orbitals, an obvious candidate 
to fill the gap is OF-DFT. The long-standing problem is 
a suitable OF approximation to the kinetic energy (KE). 
Background about the problem and a detailed descrip- 
tion of our first OF-KE functionals were reported in a 
paper addressed to the computational materials science 
community^, with a more didactic survey in Ref. i8|. The 
present analysis focuses on identifying the causes of lim- 
itations of those functionals and ways to eliminate those 
limitations. 



II. BACKGROUND SUMMARY 

Construction of an accurate, explicit total electronic 
kinetic energy density functional T[n\ — (\I>|r|^') for a 
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many-electron system in state \^) with electron num- 
ber density n is an unresolved task^ i^°'^-^ . The Coulomb 
virial theorem suggests that the task is equivalent to seek- 
ing the total energy functional itself. The Kohn-Sham 
KE is thus a more attractive target for multiple rea- 
sons. Of course, the appeal of OF-DFT predates modern 
DFT, as witness the Thomas-Fermi-Diraoi2ii^ and von 
Weizsacker^ models. There has been considerable activ- 
ity more recently. A review with extensive references is 
given in Ref. [l^- Other relevant work is that of Carter 
and co-workers; a helpful review with many references is 
Ref. ITgI. More recent developments include, for example, 
Refsr[l3ilIillHllI13 as well as our own work already 
cited. 

Distinct from most other recent efforts, our approach is 
to construct one-point, i.e., local approximate KE func- 
tional specifically for MD computations. We insist on 
constraint-based forms and parameters, that is, satisfac- 
tion of known exact results for positivity, scaling, and the 
like. We are willing, as needed, to simplify the search by 
requiring only that the functional give adequate inter- 
atomic forces, not total energies (and certainly not KS 
band structures nor general linear response). 

To summarize basics and set notation, we first note 
that except as indicated otherwise we use Hartree atomic 
units. The Kohn-Sham2^ kinetic energy Tg, the major 
contribution to T, is defined in terms of the KS orbitals: 

N 

= y"iorb(r)d3r. (1) 

The remainder, T — Tg, is included in the exchange- 
correlation (XC) functional E^dn]. Since successful E^c 
approximations assume this KS KE decomposition, we 
focus on Tg. This approach also evades the formidable 
task associated with the full T[n] just mentioned. 

For Ts[n] an explicit functional of n, the DFT total 
energy functional is orbital-free: 
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OF-DFT 



[n] = n[n] + E^.ln] + Eii[n] 



E^c[n]+E^ 
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(2) 



Here i?No['^] is the nuclear-electron interaction energy 
functional, -EhM is the Hartree functional (classical 
electron-electron repulsion) , and i?NN is the inter- nuclear 
repulsion. Then the variational principle gives the single 
Euler equation 



Sn{r) 



+ WKs(N;r) = M> 



(3) 



where /i is the Lagrange multiplier for density nor- 
malization J n{r)(Pr = N a.t the nuclear configuration 
Ri, R2, . . ., and vks = '5(-Enc + Eii + E^c)/Sn. The force 



on nucleus / at R/ is simply 



= -V 



ST, 



6n(r) 



VR,n(r)d3r. (4) 



The third term in Eq. ([4]) shows that the biggest error 
in the calculated force will come from the gradient of the 
approximate Tg [n] functional, because the kinetic energy 
is an order of magnitude larger then the magnitude of 
_Exc (which also must be approximated in practice). As 
will be discussed, when the development of approximate 
KE functionals focuses on forces, it is convenient to use 
Eq. ((H) with number density ?T.(r) from a conventional KS 
calculation (with a specific approximate E^c) as input. 

In constructing approximate functionals it is quite 
common to begin with the Thomas- Fermi functiona l^^'^'^ , 
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By itself, the TF functional is not an acceptable approx- 
imation, because of, for example, the Teller non-binding 
theorem^^. A more productive route for our purposes is 
to decompose Tg [n] into the von Weizsacker energy Ti 



plus a non-negative remainder, known as the Pauli term 

25.26.27.28 



T,[n]=Tw[n]+Te[n], Tg[n] >0 



(6) 



with 



|V^(r)P 
n(r) 



tvj[n{r)]cfr . (7) 



Previously we have shown^ that the non-negativity of Tg 
and tg{r), defined by 



Te = J tg{r)d^r 



te 



tn 



(8) 



is a crucial but not sufficient condition for determinin- 
ing a realistic OF-KE approximation. Details of some 
remaining issues follow. 



III. GGA-TYPE KE FUNCTIONALS AND 
THEIR LIMITATIONS 

A. Basic structure 

Pursuit of local approximations for torb(r) = 
^w[f^(I■), Vn(r)] -I- tg[n{r), Vn(r), . . .] stimulates consid- 
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eration of a counterpart to the generalized gradient ap- 
proximation (GGA) for XG22i, namely 



nGGAr 



Here s is a dimensionless reduced density gradient 



|V"I 
2nkp 



2^^1/3 



(9) 



(10) 



Ft is a kinetic energy enhancement factor which goes to 
unity for uniform density. Equation ^ is motivated in 
part by the conjointness conjecture'^", which posits that 
Ft{s) cx Fx{s) where Fx is the enhancement factor in 
GGA exchange. We showed previously that this rela- 
tionship cannot hold strictly-^, but the form is suggestive 
and useful. 

For connection with > 0, we re-express Tvv in a 
form parallel with Eq. dH). From Eqs. ^ and ITOl) . 



rw[n] = co / n5/3(r)|s2(r)d3r 



(11) 



Then Eq. ([6|) gives 

Tf^^M = TwM +coJ n'^Hr)Fe{s{r))d'r , (12) 
where 

Fgis)=Ftis)-^s^. (13) 

The final term of Eq. thus is a formal representation 
of the GGA Pauh term T^'^^. Note that the form of 
Eq. I|12p autornatically preserves proper uniform scaling 
of (see Ref.^): 



n^(r) = 7'^n(7r) . 



(14) 



Constraints that must be satisfied by the enhancement 
factors associated with any satisfactory GGA KE func- 
tional include 

t,(H;r) >0, (15) 

as well as^^i^i^ 

V0{[n];r) ^ STg[n]/Sn{r)>0, Vr. (16) 

The quantity ve is known as the Pauli potential. Con- 
straint Eq. implies the non-negativity of the GGA 
enhancement factor, Fg{s{r)) > 0. 

For a slowly varying density that is not itself small, we 
have s « 0, and it is appropriate to write as a gradient 
expansionSl: 



Truncation at second order in s gives the second-order 
gradient approximation (SGA), with the SGA enhance- 
ment factor— 
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(19) 



These forms should be exhibited by the exact functional 
in the limit of small density variation. (Though there are 
s oo constraints^, we have not used them so far.) 
For GGA functionals, vg, Eq. HH), can be written^^i^^ 

as 



,5TgGGA ^ dtg[n{r),Vn{r)] 
Sn{r) dn{r) 

After some tedium, one finds 



dtg[nir),Vn{r)] 
a(Vn(r)) ■ 
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ds'^ 



Vs • 



9Vr 



(21 



(The last line was omitted in Eq. (34) of Ref. but 
included in the actual numerical work.) A somewhat 
cleaner expression that also makes it easier to understand 
the extension we present below comes from shifting to 
the variable and defining both the reduced Laplacian 
density p 
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{2kF)^n 4(37r2)2/3„5/3 



(22) 



and one of the various possible dimensionless fourth- 
order derivatives q 

_ Vn • (VVn) • Vn Vn • (VVn) • Vn 



{2kpYn 



16(37r2)4/3„13/3 



(Note that our §2, p are denoted as p, q respectively in 
Ref. M-) Then Eq. jUI) becomes 



vf'^^s') = con^ 



Fgis') 



'-s^ + 2p 



dFg 

a(s2) 
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d^Fg 
9(s2)2 



(24) 



See Appendix A for details. 



B. Singularities 

Near a point nucleus of charge Z at the origin, the 
number density behaves to first order in r as 



Ts[n] = TtfM + hTw[n] + higher order terms. (17) 



n(r)~(l-2Z|r|) + 0(|rn 



(25) 
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as required by Kato's cusp condition^i^i^i^i^. SufB- 
ciently close to a nucleus therefore, n{r) behaves as a 
Hydrogen-like Is-electron density 



nnir) ^ exp(-2Z|r|) 



(26) 



That is, the variation of n'{r)/n{r) is equal for these 
densities sufficiently close to the nucleus, hence the form 
in Eq. ((26|l is a reasonable near-nucleus approximation^. 
Consequences of differences in the higher-order terms in 
the respective Taylor expansions of actual and hydrogenic 
Is densities are discussed in Section IV-A. Also see Ref. 
I43I for a related discussion. 

For n{r) of the form of Eq. (1^^ near r = 0, s and 
q remain finite while p — > — 4Z/(2fci?)^r. In this case, 
Eq. (gl) becomes 



„GGA 



(r^O) 



5r d{s'^) 



+ nonsingular terms. (27) 



If is sufficiently small that it is a good approximation 



to write Fp'^^ « 1 



form of rf'^^), Eq. ^ simphfics to 



(note that this is exactly the 



GGA(r->0) 



3aZ 
5r 



nonsingular terms. 



(28) 



in which case w^^^ tends to infinity at the nuclei with the 
same sign as the GGA parameter a. The small values of 

at the nuclei make this a general phenomenon. (Near 
typical nuclei (r 0), numerical experience shows that 

« 0.15, so that the small-s^ behavior of any Fe(s^) of 
GGA form is relevant there.) 

Equation (|28|l shows that purely GGA Pauli poten- 
tials have singularities in the vicinity of nuclear sites. In 
contrast, calculations using KS quantities as inputs show 
that the exact Pauli potential is finite at the nuclei (see, 
for example, Ref. [2^ as well as Fig. [2] below). Moreover, 
the positivity requirement for vg will certainly be vio- 
lated near the nuclei both for Fg'^^ and for any GGA 
form with a < 0. 

C. Positivity: tests and enforcement 

To explore these positivity constraints, we tested six 
published KE functionalsi^^i^iilii^^ that either are 
strictly conjoint or are based closely on conjointness. The 
test used the diatomic molecule SiO, an important ref- 
erence species for us. With LDA XC, we did a conven- 
tional, orbital-dependent, KS calculation as a function of 
bond length (details are in Appendix B). At each bond 
length, the converged KS density was used as input to 
the orbital- free E[n\ corresponding to one of the six Ts[n\ 
approximations. None predicted a stable SiO molecule. 
All six produced a < in Eq. hence all six have 

non-trivial violations of vg positivity, with all the effec- 
tive enhancement factors very close to that of the SGA, 
Fg{s) = 1 - 40/27s2. Details are in Refs. audi. Be- 
cause of the constraint violation, conjointness thus can, 
at most, be a guide. 



We enforced positivity of v^*^^ by particular parame- 
terization of Ft(s) forms based, in part, on the Perdew, 
Burke, and Ernzerhof (PBE)i^ GGA XC form: 



ly-l 
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PBE 



1 +Ois2 



2,3,4 
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iGXp4/ \ 



(s) = Ciil-e-"'" ) + 02(1-6-''^' ). (29) 



Our PBE2 form is the same as that used (with differ- 
ent parameters) by Tran and Wesolowski^ while PBE3 
corresponds to the form introduced by Adamo and 
Barone^®, but, again, with different parameters. Quite 
similar forms also were explored by King and Handy^ 
in the context of directly fitting a KS kinetic potential 
Vs = STg/Sn to conventional KS eigenvalues and orbitals; 
see Eq. ^ below. 

We fitted the parameters in the enhancement fac- 
tors, Eqs. to match the conventional KS inter- 
nuclear forces for various nuclear configurations of a 
three-molecule training set: SiO, H4Si04, and 11581207. 
With conventional KS densities as input, we found semi- 
quantitative agreement with the conventional KS calcu- 
lations for single bond stretching in 1148104, 11581207, 
H2O, CO, and N2. All had energy minima within 5 to 
20% of the conventional KS equilibrium bond-length val- 
ues. The latter two molecules provide particular encour- 
agement, since no data on C or N was included in the 
parameterization. Details, including parameter values, 
are m Ref. 0. 



D. Analysis of fitted functional behavior 

Despite this progress, there is a problem. Although 
the PBEz/ and exp4 forms can give Pauli potentials that 
are everywhere positive, yielding a > in Eq. (|^5|) . they 
are singular at the nuclei, in contrast to the negative 
singularities of previously published forms. For clar- 
ity about the developments which follow, observe that 
these nuclear-site singularities occur in vg, hence are dis- 
tinct from the intrinsic nuclear-site singularity of the von 
Weizsacker potential, which by operation on Eq. ([7]) can 
be shown to be 



_ STw[n] 1 
Sn{r) ~ 8 



iVnl' 



2V2n"l 



(30) 



The intrinsic singularity near a nucleus follows from Eq. 
and has the form 



Vw 



2Z 

T 



(31) 



Insight regarding the behavior of our modified conjoint 
functionals can be gained from consideration of the en- 
ergy density dT^^^^i^s^jds as a function of s for various 
functionals indicated by the generic superscript "appx" . 
This quantity comes from differentiation of the integrated 
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dT, (s)/ds 
dT, (s)/ds 
dT™'^(s)/ds - 
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FIG. 1: Energy density contributions to the Pauli term To as a 
function of s, presented as values of dTg{a)/da from Eq. (I32|l : 
shown are conventional Kohn-Sham dT^^{s)/ds (the refer- 
ence), our PBE2 functional, and the older PBE-TW GGA 
functional. Data are for the SiO diatomic molecule at bond 
length 1.926 A and are based on the density from fully nu- 
merical KS-LDA computations. 



FIG. 2; Conventional (reference) KS values for electronic 
density (scaled by the factor 47r(|«| — R/2)^, with R the in- 
ternuclear distance), Pauli term tg, Pauli potential vg, and 
enhancement factor Fg, calculated for points on the internu- 
clear axis using KS LDA fully numerical orbitals for the SiO 
molecule; Si at (0, 0, -0.963)A, O at (0, 0, -F0.963)A. 



contribution Tg^^^{s) of the region s(r) < s to the kinetic 
energy: 

J s(r)<s 

= j'^ds J t'^'P'"'{[n];r)6{s ~ s{r))dh . (32) 

Figure [1] shows dTg^^^{s)/ds for the SiO molecule at 
bond length R — 1.926 A (slightly stretched). Values 
are shown for our recent parameterization PBE2 that re- 
spects positivity, the Tran-Wesolowski parameterization 
of the same form (PBE-TW), and for the exact, orbital- 
dependent KS Pauli term calculated froro^^ 

= ^o.b-Q^-iv^n) , (33) 

where tg and torb are defined in Eqs. ([U and ^ respec- 
tively. Recall that the exact value of tg must be non- 
negative^^. For clarity, note also that while torb can be 
negative, the equivalent form 

ts = iorb + ^V^?! (34) 

is positive definite^ ^. 

Figure [T] also shows that both approximate functionals 
closely resemble the exact KS kernel for 0.24 < s < 0.38. 
But both of them have a much larger second peak around 
s « 0.5. In contrast, the exact functional actually has a 
long low region before a second peak at s » 0.9. Our 
PBE2 approximate functional mimics the true second 



peak via a too-strong third peak while the conventional 
GGA PBE-TW functional has a spurious minimum at 
this point. Moreover, the PBE-TW Pauli term goes neg- 
ative for all s > 0.82. In addition, we see from Fig. [1] 
that the KS kinetic energy is nearly totally determined 
by the behavior of Fg over a relatively small range of s, 
approximately 0.26 < s < 1.30 for the SiO diatomic. The 
asymptotic regions (s — > and s — > oo) do not play a sig- 
nificant role. The range 0.26 < s < 0.9 has the highest 
weight of contribution (the highest differential contribu- 
tion). As an aside, we remark that PBE2 overestimates 
the KE presumably because it was fitted purely to forces 
without regard to total energies. 

Figuresdl [21 and|3]provide comparisons of the reference 
tf^, vf^, and F^^ (where F^^ = tf^/corV''^) with the 
corresponding quantities for the PBE-TW and PBE2 ap- 
proximations. The values are along the internuclear axis 
of the SiO molecule with internuclear separation 1.926 
A. The KS Pauli potential was calculated using the ex- 
act orbital-dependent expression 

where 4>i and e.i are the occupied KS orbitals and eigen- 
values respectively. Equation ([35]) is obtained in a way 
similar to that used in Ref. [2^; see Ref. 

In Fig. H all three KS quantities, tf^, vf^, and F^ 
are everywhere non-negative, as they must be. Observe 
that is finite at the nuclei and has local maxima in 
positions close to the inter-shell minima of the electronic 
density. 
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FIG. 3; As in Fig. [2] for the PBE-TW conjoint approxi- 
mation. Lower left panel: region near Si-center; Lower right 
panel: region near 0-center. 
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FIG. 4: As in Fig. [2]for the PBE2 modified-conjoint approx- 
imation. 



so the PBE-TW misbehavior might be resolved by such 
an addition. See discussion below. Figure [3] also shows 
that vg for the PBE-TW functional has very sharp neg- 
ative peaks exactly at the nuclear positions, in accord 
with Eq. ([28)1 . 

Figure [4] shows that the modified conjoint PBE2 ve re- 
spects the positivity constraint everywhere. Inter-nuclear 
forces in the attractive region are described at least qual- 
itatively correctly as a result. The PBE2 potential is still 
divergent at the nuclei in accordance with Eq. ([28]) . 

As an aside we remark on two small computational is- 
sues. First, the absence of nuclear-site singularities in the 
computed correct Pauli potentials v^^ might be argued 
to occur because the computed density does not have pre- 
cisely the proper nuclear site cusp, i.e., does not strictly 
obey Eq. (|25p . However, our numerical results are consis- 
tent with those in Ref. [2^. Those authors used numerical 
orbitals^"'^ which presumably satisfied the cusp condition 
approximately. More importantly, if a specific numeri- 
cal technique gives a KS density and associated KS ts 
that produce a non-singular vg, then an approximate vg 
evaluated with the same density should not introduce 
singularities. Second, the reader may notice that Fig. 2] 
shows small negative values for tg far from the bonding 
region of the molecule. This behavior is due to numerical 
imprecision associated with computation for extremely 
small values of the density, and makes no appreciable 
contribution to the kinetic energy. 

Finally, an insight to the harm of excess positivity of vg 
can be seen by examining the dependence of Tg [n] upon 
Vg . From a known virial relation^^ we have 



T9[n]^^ J vgi[n];r){3 + r-V)n{r)d^r. 



(36) 



Any spurious singularities of Vg^"^^ at the nuclei clearly 
will cause special problems in overweighting the inte- 
grand. 



IV. BEYOND GGA-TYPE FUNCTIONALS 

The preceding analysis makes clear the need for more 
flexible functionals than the forms of Eq. (I^^)) . In par- 
ticular, nuclear site divergences of vg are unavoidable 
for all purely GGA-type functionals e.g., GGA, GGA- 
conjoint, and modified conjoint KE functionals; recall 
Eqs. (P7)) and (pS)) . Additional variables and constraints 
upon them are required to eliminate the singularities. 



A. Reduced derivatives of the density 



In contrast, the energy density of the PBE-TW Pauli 
term and corresponding enhancement factor have nega- 
tive peaks in the inter-shell regions, violations of the non- 
negativity constraint for tg. However, addition to tg^^^^ of 
any multiple of a Laplacian term V^n would change only 
the local behavior without altering the value of Tg^^^[n], 



Consider again the gradient expansion of Ts[n], Eq. 
(HZ]) (see Refs. m, mUll for details), which we recast 
as 

T,[n]= J {io(N;r)+i2(N;r)+<4(M;r) + ...}d3r. 

(37) 
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Here to is as in Eq. t2 = (1/9) ^w, and 



t4{[n];r) 



,5/3 



(r) 



/V2n(r)\2 
Vn5/3M/ 



/ Vn(r) \2/V n(r)\ 1 / Vn(r) \ 
8U4/3(r)y U5/3(r)/ + 3Vn4/3fr)/ 



540(37r2)2/3 
Vn(r) \2/V2n(r) 



n5/3(r) 
3' " 



/3(r). 



(38) 



The sixth-order term, dependent upon n, |Vn|, V^n, 
|VV^7i|, and V^rt, is given in Ref . [52. 

As is "well-known, in a finite system (e.g. molecule), 
a Laplacian-dependent term V^n affects only the local 
behavior of the kinetic energy density. Arguments have 
been advanced for and against including such Laplacian 
dependence in the KE density functional; for example, 
see Refs. 113, [13, [H. Recently Perdew and Constantii>i^ 
presented a KE functional that depends on V^n via a 
modified fourth-order gradient expansion. Though not 
stated that way, their functional obeys the decomposition 
of Eq. It is intended to be universal (or at least very 
broadly applicable), whereas we are focused on simpler 
functionals that require parameterization to families of 
systems. The Perdew-Constantin form involves a rather 
complicated functional interpolation between the gradi- 
ent expansion and the von Weizsacker functional. They 
did not discuss the corresponding potential vg nor Born- 
Oppenheimer forces. And they characterized the per- 
formance of their functional for the energetics of small 
molecule dissociation as "still not accurate enough for 
chemical applications" . So we proceed rather differently. 

Rearrange the foregoing gradient expansion into Tw + 
Te form (recaU Eq. I®): 



Te[n] = 



con5/3(r) (^l-^s^^ + h+te + 
ta(l-^sA + ta-^s^ + t4 + te -+ 
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t'-^\[n];r)+t);'i[n];r)+t);>i[n];r) 



.(2) 



(4), 



where 



4°)([n];r) = to(N;r) 



(39) 
(40) 



where s, p are as in Eqs. ([10]) and (|22p respectively. The 
first two terms of the expansion yield the SGA enhance- 
ment factor already discussed 



^SGA ^ ^(0) 



a2S 



(44) 



with a2 — —40/27. The fourth-order (in highest power 
of s] term is 



( 4-1 AO o 

Fa = 045 -I- b2P + C21S p , 



(45) 



with coefficients a4 = 8/243, 62 = 8/81, and C21 — —1/9. 

Rather than retain those values of 02, 0,4, 62, C21, we 
instead treat them as parameters and seek values or rela- 
tionships among them which would yield a non-singular 
vg through a given order. (Corresponding improvement 
of Thomas-Fermi theory by imposition of the nuclear 
cusp condition was introduced in Ref. 56.) 

Functional differentiation of each term in Eq. (|39p gives 



the formal gradient expansion vg = Vg^"' ~^^g^"^ + ' ' ' , 
where Fg'^^^ (shown below with its arguments suppressed 
for clarity) is a function of s^, p, and in principle, higher 
derivatives of n(r): 



„(2») 



(r) = io(N;r) 



3n(r) 



-F, 



(?(s2) dn{r) 



dp 

dp dn{r) 



V- to(N;r) 



9(s2) dVn{r) 



to([n];r 



dF, 



dp 



dp 9V^n(r) 



.(46) 



The ellipses in Eq. (|46)) correspond to additional terms 
that are needed only if Fg'^^'^ depends upon derivatives 
other than s and p. 

After manipulation (see Appendix A), one obtains the 
potentials corresponding to the enhancement factors in 
Eqs. dm) and dlil): 



Con 



2/3 



a2S 



2a2P 



(47) 



and 



if (H;r)=io(N;r) 
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4^)(M;r)=to(N;r)[^(: 



„2 9 2, 



, n s jD H — s 

81V 8^3 



(41) 



(42) 



Each term, Eq. ((40])-([42|), of Eq. ((39j) can be put straight- 
forwardly into a GGA-like form: 

tf)(M;r) = io(N;r)i^fH5,p,...), (43) 



,(4) 



CqU 



2/3 



1104 + —C21 I s' 



80 



(562 + 2c2i) - ( 4a4 - y 62 ) s^P 



804 + y C21 ) g 



20 



b2q' + 2h2q" + 2c2iq" 



(48) 
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Here q is as in Eq. p3|) and q' , q" , and q'" are other di- 
mensionless fourth-order reduced density derivatives de- 
fined as 



Vn • VV^n 



(2fcF)4n2 16(37r2)4/3„io/3 



(2A:p')4n f 6(37r2)4/3„7/3 ' 
VVn : VVn VVn : VVn 



(2/cF)4n2 16(37r2)4/3„10/3 



(49) 
(50) 
(51) 



The operation denoted by the colon in the numerators of 
q'" isA:B^j:.^^A,B,,. 

At Eq. (|77|) we have aheady pointed out that an 
enhancement factor of SGA form, specifically, that of 
Eq. (l44l) . produces a Pauli potential 



,(2), 



3 Za2 
5 r 



-nonsingular terms. 

(52) 

This exhibits the 1/r SGA Pauli potential nuclear sin- 
gularity already discussed; we return to this point in a 
moment. 

For the fourth-order enhancement factor, Eq. (|45p . we 
again note that s and q are non-singular near the nucleus, 
while 



lims2(r) = ZV[37r2n(0)]2/3 

r— +0 



(53) 



With a density of the form of Eq. ((26|) , Eq. (|48)l thus gives 
the near-nucleus behavior of the fourth-order potential as 



Co 



16[97r47i(r)]2/3 
32^3 



^(562 + 3c2i) 
3r^ 



9r 



(18a4 -I- I7b2 + 18c2i) 



nonsingular terms. 



(54) 



The singularities in l/r2 and 1 /r can be removed by re- 
quiring that the numerators of the first two terms of Eq. 
(f54| both vanish, or equivalently 



C21 



3 



^ 18 



(55) 



In the spirit of the GGA, we are led to defining a 
fourth-order reduced density derivative (RDD) as 



which yields a Pauli potential with finite values at point 
nuclei. Clearly it is not the only K4-dependent enhance- 
ment factor with that property. So, we seek Fg^\K4) 
functional forms which are more general than Eq. (j57p 
and which give a positive-definite, non-singular vg. 

At this point, it is prudent to consider how many terms 
in the Taylor series expansion of the density Eq. (|26p are 
relevant for the cancellation of singularities in Eq. (I54p . 
The answer is four terms: n{r) c>c 1 — 2Zr + 2Z^r^ — 
{A/3)Z^r^ . That is, the singularities will reappear for a 
density of the form of Eq. ([^5]) if the second- and third- 
order terms differ from those defined by a Hydrogen-like 
density expansion, e.g., Eq. ([26]) . Thus, the foregoing 
cancellation fails for a density with power series expan- 
sion n(r) cx l~2Zr-{4/3)Z^r^ + . . .. This fact wiU limit 
applications of simple K4-based KE functionals to those 
densities which have precisely Hydrogen-like behavior up 
to fourth order. 

It is necessary, therefore, to consider other candidates 
for RDD variables which would provide cancellation of 
singularities for the density Eq. (|25p independently of hy- 
drogenic higher-order terms in the Taylor series expan- 
sion of the density. The observation that K4 ^ 0(V'*) 
suggests that the effective or operational order of V in 
such a candidate variable should be reduced to second 
order. This in turn suggests a candidate variable, still 
based on the fourth-order gradient expansion Eq. (|42p . 
namely 



(4-2) 



62^2 



C2is2p, 



(58) 



(compare Eq. (|45| ) . Now consider a density of the form of 
Eq. (|25p but with arbitrary first- and higher-order near- 
nucleus expansion coefhcients. 



n{r) ^ {1 + Cir + C2r^+C3r^). 



(59) 



Following the same lines as those used to reach Eq. (|54p . 
one finds 

Vg'^ '^\r) ^ ~i= — 1^ nonsingular terms . (60) 
V ^2 ^ 

The singular term would be eliminated by the choice 
C21 = 0. The cancellation is universal in that it does not 
depend on the density expansion coefhcients, Ci (while 
the singular term prefactor and non-singular terms do, of 
course, depend on those expansion coefficients) . Hence a 
candidate RDD variable (denoted as K4) which provides 
cancellation of singular terms in the Pauli potential could 
be defined as 



K4 



52p2 



62 > 0. 



(61) 



4 18 2 30 2 frr.. 

Ki = s " P- (56) 

This RDD with Eq. (I45p gives an enhancement factor 

Flf'\K.4) ^ a^Ki , (57) 



Note that this form is manifestly positive. 

This RDD can be used to construct a variety of en- 
hancement factors to replace Eq. flS]) for the fourth- 
order approximation to the Pauli term, for example 
^9(^4) — a^iii. This simplest enhancement factor cor- 
responds to a Pauli potential with finite values at point 
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nuclei but clearly it is not the only K4-dependent one with 
that property. Any linear combination of non-singular 
enhancement factors (including the simple Fg ^ I) also 
will be non-singular. A combination of two PBE-like 
forms (see Eq. (|68| below) is also non-singular, as can 
be checked analytically for any density with near-nucleus 
behavior defined by Eq. ([59]) . hence also Eqs. ((25l) . (f26l) . 

There are, of course, K4-dependent functionals that 
yield a divergent potential, e.g. 



A.4 



(62) 



so one must be cautious. 

Regarding the second-order forms, Eq. ((52|l shows that, 
short of complete removal of the term from Fg^^, we 
cannot cure the singularity in Vg'^^. There is no direct 

analogy to the removal of singularities in Vg^"^ just dis- 
cussed. Instead, in parallel with Eq. (|56p or Eq. ((6T|) . we 
introduce a second-order RDD 



K2 = S 



biP, 



(63) 



with 61 to be determined. Then, in analogy with a PBE- 
type enhancement factor, we can define an enhancement 
factor dependent only on second-order variables as 



K2 



1 + aK2 



(64) 



For it, the near-nucleus (small r) behavior of the Pauli 
potential is 



61 



0(r) 



(65) 



(2) 

with constants > which depend on the specific 

density behavior being handled. 

The RDDs considered thus far are combinations of 
powers of s and p which ensure cancellation of nuclear 
cusp divergences in vg. Thus, we define a class of approx- 
imate KE functionals, the reduced derivative approxima- 
tion (RDA) functionals, as those with enhancement fac- 
tors depending on the RDDs 



io(M;r)Fe(K2(r),K4(r)) (fr, 



(66) 

{to[n] from Eq. ([5])) with non-divergent Pauli potentials 
as a consequence of constraints imposed on the coeffi- 
cients in the RDDs. This route of development of KE 
functionals is under active investigation; see below. 

For insight. Figure [5] shows the behavior of the RDD 
K4 along the SiO internuclear axis for four values of 62- 
The behavior in the vicinity of the Si atom is shown. 
Both the s and p variables have four maxima which lie 
close to the intershell minima in the density. Increasing 
the value of 62 increases the height of the corresponding 
maxima for K4 RDD (because the contributions from the 
p-maxima increase). 



1 ' 

4ii{iz|-R/2)^n/10 
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' 1 
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FIG. 5: The fourth order, ^4 reduced density derivative for 
different values of 62 along the internuclear axis z for the SiO 
diatomic molecule near the Si atom: Si at (0,0, -0.963) A, O 
at (0, 0, -^0.963) A. Variables s and p, are shown for compari- 
son. 



One of the peculiarities is that the reduced den- 
sity Laplacian p is divergent at the nucleus and, as 
a consequence, K4 itself also is divergent, even though 
it generates a non-divergent vg. This divergence will 
not affect the KE enhancement factors provided that 
liniK^^oo ^e('«4) = Constant. One of the advantages of 
the K4 variable is its positiveness everywhere (by defini- 
tion). 



B. Parameterization of a RDA functional 

In addition to the positive singularities, another limita- 
tion of our earlier modified-conjoint type KE functionals 
was the inability to parameterize them to provide both 
forces and total energies simultaneously^. Given the em- 
phasis on MD simulations, parameterization to the forces 
was the priority. With the spurious repulsive singulari- 
ties removed from RDD functionals, the question arises 
whether total energy parameterization can be used and, 
if so, if it is beneficial. The usual energy fitting criterion 
is to minimize 



uje 



E 



KS 



OF-DFTI 



(67) 



over systems (e.g., atoms, molecules) and configurations 
(e.g., diatomic molecule bond length) indexed generically 
here by i. When the parameter adjustment is done for 
fixed-density inputs (i.e., conventional KS densities as 
inputs), this total energy optimization is equivalent to 
optimization of the Tg functional. We did this for deter- 
mination of the empirical parameters for the new RDA- 
typc functionals F^^^ = Fg{ki). 

Since Fg ~ I (or any constant in general) also yields a 
non-singular Pauli potential, we can form K4-dependent 
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enhancement factors which resemble GGA forms and 
thereby enable connection with the modified conjoint 
GGA functionals discussed already. One form which we 
have begun exploring (see below) is 



-A. 



K4 



1 + /3iK4 
j 



K4 



1 + (32^4 



(68) 



Ai and (3i are parameters to be determined. Even this 
simple form has two desirable properties: (i) the corre- 
sponding vg is finite for densities with the near-nucleus 
behavior defined by Eq. ([55]) . hence also Eqs. (^5)1 or 
((26)) (this has been checked by explicit analytical calcu- 
lation); (ii) the divergence of K4 near the nucleus (see 
Figure E]) cancels in Eq. ([68]) (lim^i^oo F^^^^^^\k4) — 

Ao + AilPl + A2lPi). Positivity of pf^^^'^^ depends 
on the parameters Ai and must be checked for any given 
determination of their values. 

After limited exploration, we used i=2, j—A. Again 
because the motivating materials problem was brittle 
fracture in the presence of water, our choice of train- 
ing sets tended to focus on SiO. We used two molecules 
with Si-Q bonds and two closed shell atoms, AI = 
{HgSi207, £[48104, Be, Ne}, with a set of six bond lengths 
for each molecule. That is, for the H6Si207 one of 
the central Si-0 bond lengths was changed, i?(Sii- 
Oi)={1.21, 1.41, 1.61, 1.91, 2.21, 2.81} A. For H4Si04, 
the deformation was in Td mode: all four Si-0 bonds 
were changed identically, i?(Si-0,)={1.237, 1.437, 1.637, 
1.937, 2.237, 2.437} A. KS-LDA densities and energies 
were the inputs (again see Appendix B for computational 
details) . Minimization of the target function defined by 
Eq. dSIl) gave 62=46.56873, Ao=0. 51775, Ai=3.01873, 
/3i=l. 30030, ^2^-0. 23118, and ^2=0.59016. A sim- 
ple check shows that the resulting enhancement factor 
^RDA(24) j.-^^ positive for all positive values of K4 (re- 
call that /{4 is positive by definition). Figures [6] and [7| 
show the p^^^^'^^^ enhancement factor as a function of 

for selected values of p and, reciprocally, as a function 
of p for selected values of s^. This is a smooth, posi- 
tive function. pf'^^'^'^^\s'^ ,p > 0.4) becames practically 
a straight line, essentially independent of s^. 

Table U displays kinetic energies for 14 molecules (four 
of them with Si-0 bonds) calculated at equilibrium ge- 
ometries by the conventional KS method and by approxi- 
mate OF-DFT functionals using the KS density as input. 
The Thakkar empirical functional'*^ was chosen as an ex- 
ample of a GGA KE functional. The Perdew-Constantin 
meta-GGA^^ "MGGA" was chosen because it, like our 
functional, is based on quantities that are at fourth or- 
der in the density gradient expansion. The results are a 
bit surprising, because the mean absolute error (MAE) 
for the RDA functional is almost a factor of two smaller 
then the MAE for the Thakkar KE, and almost ten times 
smaller then the MAE for the MGGA. Given the param- 
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FIG. 6: The RDA(24) enhancement factor as a function of 
for selected values of p. 




FIG. 7: The RDA(24) enhancement factor as a function of p 
for selected values of s^. 



eterization to a small training set containing only two 
molecules with Si-O bonds and two closed shell atoms, 
we had not expected to obtain such good transferability 
to other systems. 

Because our objective is a KE functional capable of 
predicting correct interatomic forces, one of the impor- 
tant aspects is the behavior in the attractive regions of 
the potential surface. Table |TT] shows energy gradients 
for the molecules in Table U calculated at the stretched 
bond length(s) for which the "exact" (i.e. reference) KS 
atractive force has maximum magnitude. One and two 
bonds were deformed in the water molecule (respectively 
denoted in the table as H20(1R) and H20(2R)), while 
SiH4 and H4Si04 were deformed in Td mode, and only 
one Si-0 bond was stretched in H4Si0 and H6Si2 07. 

The forces were calculated by a three-point centered 
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-74.5 



TABLE I: KS kinetic energy Ts values (in Hartrees) for se- 



lected molecules and differences (T^ 



OF-DFT 



calculated 



using a GGA (Thakkar), MGGA, and RDA(24) explicit semi- 
local approximate functionals. LDA-KS densities for LDA 
equilibrium geometries (calculated as described in Appendix 
B) were used as input. 





KS 


Thakkar 


MGGA 


RDA(24) 


H2 


1.080 


-0.022 


0.103 


0.006 


LiH 


7.784 


0.021 


0.296 


0.063 


H2O 


75.502 


-0.285 


0.318 


-0.128 


HF 


99.390 


-0.353 


0.329 


-0.148 


N2 


108.062 


-0.340 


0.300 


-0.041 


LiF 


106.183 


-0.261 


0.566 


0.086 


CO 


111.832 


-0.333 


0.300 


-0.074 


BF 


123.117 


-0.273 


0.456 


0.077 


NaF 


260.097 


-0.348 


1.295 


0.648 


SiH4 


290.282 


0.084 


3.112 


0.381 


SiO 


362.441 


-0.262 


2.825 


0.293 


H4SiO 


364.672 


-0.163 


3.338 


0.293 


H4Si04 


587.801 


-0.860 


4.133 


-0.034 


HsSiaOr 


1100.227 


-1.408 


7.968 


0.086 


MAE" 




0.358 


1.810 


0.168 



KS 

Thakkar 

MGGA 

RDA(24) 




" MAE=mean absolute error 



FIG. 8: Total energy as a function of the O-Hi distance for 
the H2O molecule (with O-H2 kept at its equilibrium value) 
obtained from a KS calculation with LDA XC and from ap- 
proximate GGA(Thakkar), MGGA(Perdew-Constantin) and 
RDA(24) functionals. The LDA KS densities (calculated as 
described in Appendix B) were used as input to the orbital- 
free functionals. 



TABLE II: Energy gradient (Hartree/A) calculated at point 
Rm corresponding to the extremum of attractive force as cal- 
culated by the KS method. Approximate OF-DFT energy 
gradients are obtained by replacing ^ by yOF-DFT^ LBA- 
KS densities for LDA equilibrium geometries (calculated as 
described in Appendix B) were used as input. 

Rm, k KS Thakkar MGGA RDA(24) 

H2 1.2671 0.164 0.029 0.112 0.005 

LiH 2.455 0.046 0.016 0.037 0.016 

H20(1R) 1.3714 0.216 -0.050 -0.073 0.249 

H20(2R) 1.3714 0.416 -0.127 -0.163 0.424 

HF 1.3334 0.232 -0.071 -0.003 0.180 

N2 1.3986 0.576 -0.349 -0.819 0.244 

LiF 2.0405 0.079 -0.019 -0.032 -0.007 

CO 1.4318 0.474 -0.248 -0.659 0.466 

BF 1.6687 0.207 -0.037 -0.118 0.204 

NaF 2.4284 0.067 -0.007 1.169 -0.008 

SiH4 1.9974 0.447 0.102 0.189 0.101 

SiO 1.9261 0.278 -0.098 -0.281 0.175 

H4SiO 2.057 0.162 -0.027 -0.086 0.151 

H4Si04 2.037 0.712 -0.278 -0.714 0.745 

H6Si207 2.010 0.194 -0.022 -0.173 0.165 



finite-difference formula. As found in our previous work 
and summarized above, the GGA functionals (with the 
Thakkar functional as the example GGA functional here) 
are generally incapable of predicting the correct sign (at- 
traction) for the force; the only molecules in Table HIl for 
which GGA predicts attraction are H2, LiH, and SiH4. 
We find the situation with the MGGA functional to be 
very similar; the predicted energy gradient has the wrong 
sign in most cases. In contrast, for all but two of the table 



entries the RDA(24) functional predicts the correct sign 
of the gradient. In many cases (H20(1R), H20(2R), HF, 
CO, BF, SiO, H4SiO, H4Si04, HgSisOr) it yields values 
very close to the reference KS results. 

Figure [5] shows the energy for the water molecule as 
a function of the 0-Hi bond length. Again, neither the 
GGA nor the MGGA curve exibits a minimum. This is a 
case for which the new RDA(24) functional behaves rel- 
atively poorly. It does reproduce a minimum, but at too 
large a bond length, while in the tail region its curve goes 
almost flat. Thus the structure predicted by RDA(24) 
would be more expanded than the correct value and the 
attractive force in the tail region would be significantly 
underestimated . 



C. Atomic analysis of the RDA(24) functional 

For analysis of the new functional, we calculated the 
Pauli potential near the nucleus (r — > 0) for the Be atom 
using a simple H-like density. A single-^ Slater orbital 
density with exponents Cis = 3.6848 and (2s — 0.9560 
taken from Ref. 57 near r = has the following Taylor 
series expansion: n{r) sa 415.0479 x (1 - 7.3971 r). The 
density n{r) — 415.0479 x exp(— 7.3971 r) has the same 
slope at r = 0. It can be used as an approximate density 
for the Be atom near r = to calculate the Pauli poten- 
tial for the RDA(24) Eq. for the GGA^^ and for 
the PBE2 modified conjoint GGA functionals. Calcula- 
tions were performed using our own Maple code. For 
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RDA(24), we find 

^RDA(24)^^ ^ 0) = 14010- 5.9025 x 10^ r 

+ 1.1271 X 10^ 7-2 + 0(r3) , (69) 
while for the GGA the result is 

prA. ^ -3.1961 

^GGA/ = + 272.95 

r 

-1319.1 r + 3257.5 + 0{r^) , (70) 
and for the PBE2 modified conjoint GGA 

^PBE2(^^0) ^ 1^ + 265.31 
r 

-1319.6 r + 3257.1 +0(r3) . (71) 

As expected, the first term in the GGA potential is di- 
vergent and negative, while the PBE2 modified conjoint 
GGA functional has a divergent but positive first term. 
The numerical coefficients for the rest of the terms are 
very close for the GGA and PBE2 modified conjoint GGA 
functionals. This closeness is consistent with the analysis 
in Ref. H. 

Figure [9] (upper panel) shows the KS-LDA Pauli po- 
tential for the Be atom. Its value at the nucleus is ap- 
proximately 3.5 Hartrees. The RDA(24) potential has 
a finite (and positive) value at the nucleus, but it is a 
strong overestimate; recall Eq. (j69p . Further compari- 
son shows that the slope of the RDA(24) Pauli poten- 
tial at the nucleus has a large negative value, whereas it 
should be very close to zero (see upper panel of Fig. [9]). 
To complete the study, three enhancement factors, KS, 
GGA and RDA(24), are shown in the same panel. Again, 
as was seen in Figure [3] for the SiO molecule, there is an 
important region where F^*^^ is negative. The RDA(24) 
enhancement factor is positive eveywhere, but is still far 
from being an accurate approximation to the KS form. 
^RDA(24) j^^^ ^ sharp minimum near r w 0.75 A which 
corresponds to the change of sign of the reduced Lapla- 
cian of the density {p) (see the lower panel) . The RDD 
K4 (shown in lower panel) also has a sharp minimum near 
this point. 

V. DISCUSSION AND CONCLUSIONS 

Success for the OF-DFT calculation of quantum forces 
in molecular dynamics requires a reliable explicit form 
for Tg. Though previously published GGA- type (conjoint 
and nearly so) KE functionals yield reasonable KE val- 
ues, they fail to bind simple molecules even with the cor- 
rect KS density as input. Therefore they produce com- 
pletely unusable interatomic forces. This poor perfor- 
mance stems from violation of the positivity requirement 
on the Pauli potential. Our first remedy was to constrain 
conjoint KE functionals to yield positive-definite Pauli 
potentials. Those functionals generate bound molecules 
and give semi-quantitative inter-atomic forces. However, 
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FIG. 9: (Upper panel): numerical KS-LDA density of the 
Be atom, KS Pauli potential and corresponding enhancement 
factor, and two approximate enhancement factors. 
(Lower panel): variables s, p, K4 (62 = 50) calculated for the 
Be atom using KS-LDA density (also shown). 



they are singular at the nuclear positions, hence severely 
over-estimate the KS kinetic energy. Examination of 
the near-nucleus behavior of both original conjoint and 
modified-conjoint GGA functionals shows that the sin- 
gularities cannot be eliminated within that simple func- 
tional form. 

Truncation of the gradient expansion, at higher orders 
in s and p, allows us to identify near-nucleus singular be- 
havior and obtain relationships among the coefhcients of 
those truncations that will eliminate such singularities. 
The resulting reduced density derivatives and related 
reduced-density-approximation functionals are promising 
for the simultaneous description of kinetic energies and 
interatomic forces. 

Two other aspects of the numerical results in Table |I] 
relate to exact constraints, hence deserve brief comment. 
First, the von Weizsacker KE is the exact Ts for two- 
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electron singlets. We have not enforced that limit, yet 
the error from RDA(24) in H2 is only 6 mHartree. Sec- 
ond, violation of iV-representability by an approximate 
rj'PP'^[n] is signaled by TJ'pp'^[7i] - Ts[n] < for at least 
one nM. Five of the RDA(24) entries in Table[T]have such 
negative differences, while MGGA has none and the ear- 
lier GGA by Thakkar has many. However, interpretation 
of those computed differences is a bit tricky, in that they 
do not correspond to the rigorous A^-representability- 
violation test but to TJ'PP''[n] - T,L°^[n]. With that 
limitation in mind, the results in Table |T] are at least 
suggestive of the notion that the MGGA and RDA(24) 
functionals are A'^-representable or, in some operational 
sense, close. It is also important to remember that the 
narrow goal is a functional that can be parameterized to 
a small training set which is relevant to the desired ma- 
terials simulations. This limitation of scope is a practical 
means for limiting the risks of non-iV-representability. In 
addition, we have many RDD forms open for exploration 
other than RDA(24). 

We close with a final word of caution. The functional 
forms we have examined to this point, e.g. Eq. (j68p . may 
be too simple to provide robust and transferable KE func- 
tionals for practical OF-DFT applications. Moreover, the 
use of RDDs as basic variables in kinetic energy enhance- 
ment factors guarantees the finiteness of the correspond- 
ing Pauli potential only for those densities which satisfy 
a generalization of Kato's cusp condition Eq. ((59|) . and 
does not guarantee the satisfaction of the non-negativity 
property, Eqs. p^ -p6 )) . The latter constraint must be 
enforced separately. Nevertheless the RDA scheme ap- 
pears quite promising and further development of it is 
underway. 



where = 4(37r2)2/3 and cq = ^{3nY^^- We also 
remind the reader that 
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recall that A : B = ^ij^ji- addition, we intro- 
duce the sixth-order reduced density derivatives 
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Finally, we write n instead of the more explicit form n(r), 
and we note that the derivatives of and p with respect 
to n, Vn, and V-n are 
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Substituting these relations into Eq. ((46|) , and restricting 
consideration to cases where Fg depends only on and 
p (for which all the terms we need to use are explicitly 
shown in that equation), we have immediately 
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At this point we remark that some terms that would 
otherwise be expected in Eq. are absent because of 
the zeros in Eq. ((75|) . 



To proceed further we need to expand the last two 
terms of Eq. (|75|) . The first of these terms expands into 



APPENDIX A. 



FUNCTIONAL DERIVATIVES 

OF Eg 



This appendix provides detail relative to derivation of 
the formulas given in Eqs. dH]), (gT]), and (gH). AU of 
them follow from an evaluation of an expression of the 

,(2») 
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generic form presented as Vg (r), Eq. That equa- Observing now that V(n ^) 



tion is a straightforward expression of the rules for the 
evaluation of a functional derivative. For clarity in what 
follows, we restate here the definitions 



8s2 

V(s2) = -^Vn 
3n 



-n ^Vn and that 
2Vn • VVn 



Vn • Vn 

^2^8/3 ' 



V^n 



to([n];r) = con5/3, (72) 



Vp 



^2„8/3 

VV^n 5 {V^n)Vn 
^2„5/3 " 3 ^2^8/3 



(77) 



(78) 



14 



the expression in Eq. ([77|) can be brought to the form 
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Continuing now to the final term of Eq. (|76p. we ex- 



pand the Laplacian, obtaining initially 
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Then, combining material from Eqs. (ffS)) . ([75]) . ([SO)) , and ([5T|) . and introducing the notations in Eqs. ([T^ - fTi]) . we 
obtain the final result, applicable for any Fg that depends only on s and p: 
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We may now specialize Eq. (|82|) to the cases needed 
in the present work. Taking first i^^^^, which has no p 



dependence, all the terms of Eq. l[82|l containing deriva- 
tives with respect to p vanish, leaving only the expression 
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previously given as Eq. p4|) . 

Turning next to the specific forms of Fg discussed in 
Section IV-A, we note tliat Fg'^^ = 1 + a2S^ is not only 
independent of p, but is also linear in s^, so dFg/d{s'^) — 
02 and d'^Fg/d{s'^f = 0. This causes v'^'^^ to have the 
form 



,,SGA 



Con 



2/3 



2p 



(83) 



which simplifies to the result given in Eq. ([47]). 

Finally, we consider Fg'^^ as given in Eq. All the 

third derivatives of Fg in Eq. ([5^ vanish; the first and 
second derivatives of Fg have simple forms. We have 
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(84) 



Equation ([M]) simplifies to the result given as Eq. 



APPENDIX B. COMPUTATIONAL METHODS 



We assess functionals by comparing results from them 
with those of conventional orbital-dependent Kohn-Sham 
calculations in the local density approximation (LDA), 
using standard methods described in, for example, Refs. 
'S^.ra.illililiiiiii. The reference molecular KS cal- 
culations were done with a triple-zeta basis with polar- 
ization functions (TZVP)^^^. All integrals were cal- 
culated by a numerical integration scheme that, follow- 
ing Becked, uses weight functions localized near each 
center to represent the multicenter integrals exactly as 
a sum of (distorted) atomic integrals. Radial integra- 
tion of the resulting single-center forms is accomplished 
by a Gauss-Legendre procedure, while integration over 
the angular variables is done with high-order quadrature 
formulas developed by Lebedev and coworkers^LZ^ with 
routines downloaded from Ref. 73. These computations 
were performed using routines developed by Salvador and 
Mayer— and included in their code fuzzy. The Vosko- 
Wilk-Nussair LDA^^ was used. 

Given the KS density, for each OF functional un- 
der study we computed the total energy £:OF-dft fj-Q^ 
Eq. (21) and the interatomic forces from Eq. The 
result is a non-self-consistent calculation which tests 
whether a given OF functional can reproduce Ts[nKs]i 
or at least y-RTs[nKs] if npcs is provided. There is no 
sense in trying to solve Eq. ([3]) with an approximate OF 
functional that cannot pass this test. 
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